function [data] = aggregateCtsTimeData_pairAdoptions(data)
	
	spell2period = data.spell2period;
	
	% Read dimensions
	K1         = data.dims.dims1{5};
	K2         = data.dims.dims1{6};
	NumSpells = length(spell2period);
	T         = double(max(spell2period));
	
	% Match each day tt to index of earliest spell of that day
	period2firstSpell = accumarray(spell2period, [1:NumSpells]', [T 1], @min); % T x 1
	% (for each day tt: this vector gives index of spell that starts at the start of day tt)
	
	%%% Simplify the data to make it at the day level (not spell level)
	data.dims.NumObs = T*K1*K2;
	data.dims.dims1{2} = T*K1;
	data.dims.dims1{3} = T*K2;
	data.dims.dims1{4} = T;
	
	% T-K1 part
	myX     = data.Xparts{2}.X;     % (NumSpells*K1) x NumX
	myX_FEs = data.Xparts{2}.X_FEs; % (NumSpells*K1) x NumX_FEs
	
	myX     = reshape(myX,     [NumSpells K1 size(myX, 2)]);     % NumSpells x K1 x NumX
	myX_FEs = reshape(myX_FEs, [NumSpells K1 size(myX_FEs, 2)]); % NumSpells x K1 x NumX_FEs
	
	myX     = myX(period2firstSpell,:,:); % T x K1 x NumX
	myX_FEs = myX_FEs(period2firstSpell,:,:); % T x K1 x NumX_FEs
	
	myX     = reshape(myX,     [T*K1 size(myX, 3)]);     % (T*K1) x NumX
	myX_FEs = reshape(myX_FEs, [T*K1 size(myX_FEs, 3)]); % (T*K1) x NumX_FEs
	
	data.Xparts{2}.X = myX;
	data.Xparts{2}.X_FEs = myX_FEs;
	clear myX;
	
	% T-K2 part
	myX     = data.Xparts{3}.X;     % (NumSpells*K2) x NumX
	myX_FEs = data.Xparts{3}.X_FEs; % (NumSpells*K2) x NumX_FEs
	
	myX     = reshape(myX,     [NumSpells K2 size(myX, 2)]);     % NumSpells x K2 x NumX
	myX_FEs = reshape(myX_FEs, [NumSpells K2 size(myX_FEs, 2)]); % NumSpells x K2 x NumX_FEs
	
	myX     = myX(period2firstSpell,:,:); % T x K2 x NumX
	myX_FEs = myX_FEs(period2firstSpell,:,:); % T x K2 x NumX_FEs
	
	myX     = reshape(myX,     [T*K2 size(myX, 3)]);     % (T*K2) x NumX
	myX_FEs = reshape(myX_FEs, [T*K2 size(myX_FEs, 3)]); % (T*K2) x NumX_FEs
	
	data.Xparts{3}.X = myX;
	data.Xparts{3}.X_FEs = myX_FEs;
	clear myX;
	
	% T part
	data.Xparts{4}.X     = data.Xparts{4}.X(period2firstSpell,:);     % T x NumX
	data.Xparts{4}.X_FEs = data.Xparts{4}.X_FEs(period2firstSpell,:); % T x NumX_FEs
	
	% Aggregate Y (take sum across spells of that day tt)
	M = sparse(T, NumSpells);
	ind = sub2ind([T NumSpells], spell2period, [1:NumSpells]');
	M(ind) = 1;
	data.Y = reshape(data.Y, [NumSpells K1*K2]);
	data.Y = M * data.Y;
	data.Y = data.Y(:); % (T*K1*K2) x 1
	
	data.logLambdaOffset_dim1 = data.logLambdaOffset_dim1(period2firstSpell,:); % T x K1
	data.logLambdaOffset_dim2 = data.logLambdaOffset_dim2(period2firstSpell,:,:); % T x 1 x K2
	
	%%% No need: logLambdaOffset_dim2 >= 0 always given that logLambdaOffset_dim2 is log(popAtRisk on provider side) in this pair adoption model.
%	tmp = repmat(data.logLambdaOffset_dim2, [1 K1 1]); % T x K1 x K2
%	tmp = reshape(tmp, [T*K1*K2 1]); % (T*K1*K2) x 1
%	data.Y(find(tmp < -1e9)) = 0;
%	clear tmp;
	
	data = rmfield(data, 'spell2period');
	data = rmfield(data, 'spellDurations');
	
	myidxes = find(data.Y(:)>0);
	data.LL_offset = -sum(gammaln(1+data.Y(myidxes)));
end
